/*********************************************************************************
 * Software License Agreement (BSD License)
 * Copyright (c) 2012, Carlos D. Correa (info@ngraph.org)
 * All rights reserved.
 * 
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 * 
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. Neither the name of netzen nor the
 *    names of its contributors may be used to endorse or promote products
 *    derived from this software without specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 *  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *********************************************************************************/
#include "timing.h"

#include "node.h"
#include "sensitivity.h"
#include <sstream>
#include <math.h>

#include <assert.h>
#include <algorithm>

Sensitivity::Sensitivity(Controller *controller):AnalysisTool("Sensitivity","Moving Least Squares",controller) {
  setPropertyFloat("Radius", 0.1);
  setPropertyStringList("Property", "Default", "node", true);
  setPropertyStringList("dX", "Default", "node", false);

}

void Sensitivity::compute() {
  Graph *graph = controller->getGraph("graph");

  timestamp t1 = now();

  string props = getPropertyString("Property");
  istringstream iss(props);
  string token;
  vector<string> properties;
  while(getline(iss, token, ':')) {
    if(token.length()>0) {
      properties.push_back(token);
    }
  }


  float radius = getPropertyFloat("Radius");
  string dX = getPropertyString("dX");

  float dprogress = 1.0/(float) properties.size();
  float progress = 0;
  for(unsigned int k = 0 ; k < properties.size(); k++) { 
    string dY = properties[k];
    string strDeriName = dY+"/"+dX;
		ElementIterator *it = graph->createNodeIterator();
		for(it->begin();it->hasMoreElements();) {
			Node *node = (Node*) it->next();
      progress+=(float) 1.0/(float) graph->numNodes()*dprogress;
      if(listener) listener->setProgress(progress);

			float x = graph->getAsFloat(node, dX);
			float y = graph->getAsFloat(node, dY);
			float xn = graph->getAsFloat(node, dX, true);
			float yn = graph->getAsFloat(node, dY, true);

      float sumdxdx = 0, sumdydx = 0;
			float sumdx = 0, sumdy = 0;

      float sumw = 0;
			int neighbors = 0;
			ElementIterator *it2 = graph->createNodeIterator();
			for(it2->begin();it2->hasMoreElements();) {
				Node *ndj = (Node*) it2->next();
				float xj = graph->getAsFloat(ndj, dX);
				float yj = graph->getAsFloat(ndj, dY);
				
				float xjn = graph->getAsFloat(ndj,dX,true);
				float yjn = graph->getAsFloat(ndj,dY,true);
				float dis = ((xn-xjn)*(xn-xjn)+(yn-yjn)*(yn-yjn));
				
				bool inRange = false;
				if(dis<radius*radius){
					inRange = true;
				} else {
					inRange = false;	
				}            
				if(inRange) {	 
					neighbors++;
					float w = exp(-5*5*dis/radius);
					float dx = xj-x; //records[j].values[xaxis] - records[i].values[xaxis];
					float dy = yj-y; //records[j].values[yaxis] - records[i].values[yaxis];
					sumdxdx+= w*dx*dx; 
					sumdydx+= w*dy*dx; 
					sumdx+=w*dx;
					sumdy+=w*dy;
					
					sumw +=w;
					
				}//endof if
			}//endof for
			graph->destroyIterator(it2);
		
			sumdxdx = (sumw>0)? sumdxdx/sumw:0;
			sumdydx = (sumw>0)? sumdydx/sumw:0;
			sumdx = (sumw>0)? sumdx/sumw:0;
			sumdy = (sumw>0)? sumdy/sumw:0;
			
			//save deri
			float fDeri = (sumdxdx>0)? sumdydx/sumdxdx: 0 ;
			
			//fDeri =  (sumw*sumdxdx - sumdx*sumdx>0)? (sumw*sumdydx - sumdx*sumdy)/(sumw*sumdxdx - sumdx*sumdx): 0;
				
			node->setPropertyFloat(strDeriName,fDeri);
    }//end iterator
		graph->destroyIterator(it);
    controller->addDataProperty("node", strDeriName,PROPERTY_FLOAT);
    graph->computeNodeMinMaxFloat(strDeriName);
  } // end properties

  if(listener) listener->done();

  timestamp t2 = now();
  printf("[SENSITIVITY] Time  %f sec\n", t2 - t1);
}


